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During this period investigations were carried out in three areas: 

1. A simulation study was conducted to investigate very early stages of a nucle- 

ation and growth mechanism taking place on diamond surfaces. The ener- 
getically most favorable binding sites were calculated along with the binding 
energies for up to three carbon adatoms deposited on (100) and (111) planes 
of a diamond crystal. Bindings for carbon atoms deposited on the (100) plane 
were found to be much stronger than the (111) case. Based on purely energetic 
considerations, results obtained in this investigation indicate that the nucle- 
ation of carbon atoms should take place more favorably on the (100) plane 
than on the (111) plane. Calculated high energy sites for adatoms on the 
(100) plane were found to be commensurate with atomic sites within the sur- 
face region of the crystal. Computational details along with a short discussion 
of the results are presented in Appendix 1. 

2. Calculations were carried out in this investigation to analyze energies and 

structural properties for (2x1) restructured patterns of the (100) surface 
of diamond. Two different model functions developed recently for carbon 
were employed in calculations. Three differently reconstructed (2x1) pat- 
terns (namely, one dimer model and two different missing row models) for 
the exposed surface layer, were taken into consideration. Results obtained in 
this investigation support the stability of the dimer model over the unrecon- 
structed (1x1) surface. Calculated energies and structural properties for the 
dimer model were found to be in agreement with values obtained in other the- 
oretical calculations. While one of the potential functions, in accordance with 
common belief, predicted the dimer model for the (2x1) restructured surface 
to be the energetically most stable structure, the other function produced one 
of the missing row models as the most favorable form. Appendix 2 includes 
details of the calculated results and a discussion on the configurational aspects 
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of (2x1) restructured surfaces. 


3. In this part, a crack propagation process for the graphitic basal plane, was 
investigated using a molecular dynamics technique. Interactions among car- 
bon atoms in the system were calculated using the Tersoff potential energy 
function. This function is based on two- and three-body interactions and pro- 
duces many macroscopic properties of carbon species correctly. In addition to 
individual atomic energies, stress tensors for every atom in the system, were 
also calculated. Simulations were performed for a system containing up to 
1600 carbon atoms and bearing an initial crack. After an equilibration period 
the system was loaded uniaxially in a direction perpendicular to the existing 
crack with small incremental elongations. Further equilibrations caused the 
crack to propagate. During this crack propagation process, atoms located at 
the tip of the crack were found to have very high stress values. This result is in 
general agreement with other reported calculations and with existing theories 
based on macroscopic considerations. However, no dislocation formation in 
front of the advancing crack was found. This outcome was attributed to the 
short range nature of the potential function employed in this study. Results 
are presented in Appendix 3. 
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Appendix 1. Modeling for Early Stages of Diamond Growth 


An atomic level understanding of the very early stages of a nucleation process 
is highly desired today in different disciplines related to crystal growth [1-3]. In 
this study a simulation calculation was performed to investigate early stages of the 
diamond growth mechanism. Binding energies and binding sites were determined for 
a number of carbon atoms deposited on (100) and (111) faces of a diamond crystal. 
Calculations were carried out for the low temperature limit and throughout this 
investigation model potential functions developed recently for carbon species were 
used. For comparison, two different model functions, based on two and three-body 
interactions, were considered here. The first potential is the Tersoff function. It has 
been shown that this function is able to reproduce correctly various bulk properties 
of diamond and of the graphitic plane [4]. The second potential used in this study 
is the Brenner function [5] which is analytically similar to the Tersoff function, but 
parametrized differently. This function also, has been shown to produce acceptable 
results for properties of bulk diamond and the basal plane of graphite, as well as 
for some properties of small carbon clusters [5,6]. 


Binding energies per adatoms, were calculated as: 

- E t = - E°] 

m 

where m is the number of adatoms, E° denotes the total equilibrated energy of the 
system of N particles with an exposed surface, and E is the total relaxed energy 
of the same system with m adatoms deposited on the surface. In this investigation 
(100) and (111) index planes of a diamond crystal were taken into consideration. 
Surfaces were first generated as abrupt terminations of bulk diamond oriented to 
produce the desired surface planes. Systems bearing these faces were then relaxed 
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using a molecular dynamics procedure. Periodic boundary conditions in two di- 
rections (parallel to the exposed surface) were imposed on the system to provide 
the continuity of the surfaces. At low temperature these planes remain in (1 x 1) 
patterns after the relaxation, but results indicate that for both surface planes the 
top interlayer spacings contract considerably, while the second interlayer spacings 
exhibit a moderate expansion [7]. Next, carbon atoms were positioned on these 
fully relaxed surfaces in various configurations and the whole system (now includ- 
ing adatoms) was reequilibrated with the same procedure. Up to three adatoms 
were deposited on these surface planes and for each case a number of binding sites 
were considered in calculations. Only the energetically most favorable sites are 
reported here. 


A top view of the (100) surface plane of diamond is depicted in Figure 1 which 
shows lattice sites for atoms located in the top four layers. Calculations using 
both Tersoff and Brenner functions produced an upper bridge site (which is the 
mid-point between two neighboring upper layer atoms) as the energetically most 
favorable site for a single adatom. One of these sites is denoted by the letter A in 
Figure 1. These upper bridge sites are located exactly above the atoms in the fouth 
layer and they represent exact positions for adatoms forming a complete monolayer 
coverage over the (100) plane of diamond. Calculated values for binding energies 
are given in Table 1. Equilibrated positions for the adatom in this case, were found 
to be about 0.9 and 0.85 A above the surface plane for the Tersoff and Brenner 
functions, respectively. For two adatoms deposited on the (100) plane the lowest 
energy configuration using the Tersoff function was found to be the two neighboring 
upper bridge sites denoted by A and D in Figure 1. For the Brenner function, on 
the other hand, the two opposing upper bridge sites indicated by A and B are 
more favorable. In this case also, after the equilibration, adatoms were found to be 
located approximately 0.9 and 0.85 A above the surface plane for the Tersoff and 
and Brenner functions, respectively. In the case of three adatoms, both functions 
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produce sites A, B and C (see Figure 1) as the energetically most favorable sites. 
The adatom at site C in this case is located about 1.75 A above the surface plane for 
the Tersoff function and about 1.65 A above the surface for the Brenner potential. 
In both cases, the deposited t rimer in its fully relaxed configuration forms an obtuse 
isosceles triangle with an apex angle of approximately 113 degrees. 

In the case of the (111) plane, lattice sites for atoms located in the top four 
layers are shown in Figure 2. For an unrelaxed surface, only one type of interplan ar 
spacing exists for the (100) plane. For the (111) plane, however, atomic layers come 
in pairs and form two different types of interlayer spacings. For an ideal case, the 
first interlayer spacing, £12 is just j of the second interlayer spacing, 82 3 , which is 
equal to the nearest neighbor distance in the crystal. In (111) plane, a considerable 
amount of shrinkage takes place in £12 upon relaxation [7]. Calculations using both 
the Tersoff and Brenner functions, produced the energetically most favorable site 
for a single adatom on the (111) surface plane at a top position. The location for 
one of these binding sites is denoted by the letter A in Figure 2. The optimized 
position of an adatom was calculated as 1.49 and 1.38 A exactly above the surface 
atom for the Tersoff and Brenner functions, respectively. In both cases, calculated 
results also indicate that the surface atom at site A moved upward for a few tenths 
of an angstrom above the surface plane. For a dimer deposited on a relaxed (111) 
surface, calculations based on Tersoff function produced adatom positions near two 
neighboring top surface atoms as the energetically most favorable sites (A and B as 
shown in Figure 2). In its fully relaxed configuration adatoms of the deposited dimer 
are about 1.62 A above the top layer. Adatoms in this case are not exactly above 
the sites A and B but displaced symmetrically toward each other. Equilibrated 
adatoms were found to be aligned parallel to the AB line with an interatomic 
distance of 1.48 A. In this case also, surface atoms at sites A and B relaxed upwards 
for approximately two tenths of an angstrom. For the Brenner function the relaxed 
configuration of a dimer is somewhat different. In this case, adatoms were found 
to be aligned above the line connecting the sites A and C. One of the adatoms was 
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located between A and D about 1.67 A above the surface plane while the other 
adatom was found to be between C and D about 1.92 A above the surface. In this 
fully equilibrated configuration the distance between two adatoms was calculated 
as 1.36 A. Surface atoms at sites A and C were found to be relaxed upwards about 
0.3 and 0.6 A above the surface plane, respectively. In this case it is worthwhile to 
indicate that the atom at site C was located in the second layer before the relaxation. 
In the case of three adatoms, both functions produced comparable results. Adatoms, 
in their energetically most stable configuration, are located near the sites A, D and 
C (see Figure 2). For the Tersoff function adatom positions were found to be about 
1.73, 2.53 and 1.99 A above the surface plane, while the Brenner function produced 
1.71, 2.43 and 1.98 A, respectively. In both cases, adatoms form an obtuse isosceles 
triangle with an apex angle of approximately 125 degrees and bond distances were 
calculated as 1.46 and 1.31 A for the Tersoff and Brenner functions, respectively. 
After relaxation, in this case also, surface atoms at sites A and C were found to be 
displaced upward approximately 0.2 and 0.5 A for the Tersoff function and about 
0.3 and 0.6 A for the Brenner function. 


Calculated results for E & reported in Table 1 indicate that binding of adatoms 
on the (100) surface is energetically more favorable than binding on the (111) plane. 
At least up to m = 3, high energy adsorption sites on the (100) face were found 
to be commensurate with the atomic sites on the surface of the crystal. However, 
high energy positions for adatoms (of dimers and trimers, in particular) deposited 
on the (111) plane, are only partially commensurate with the lattice sites. Both 
potential functions produced comparable results. Closer analysis of structures of 
equilibrated trimers, reveals that at its very early stages a nucleus formed on the 
(100) plane is more likely to have a diamond character than one forming on the 
(111) surface. Based on purely energetics considerations, results obtained in this 
study indicate that carbon atoms nucleate on the (100) surface more readily than 
on the (111) plane. In comparing these results with experimental findings, however, 
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extreme care should be exercized. Calculations carried out here employing model 
functions reflect energetics of the nucleation process for systems containing only 
carbon atoms. No consideration was given in this investigation to the entropic or 
kinetic aspects. 
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Table 1 . Calculated binding energies for carbon atoms deposited on (100) and 
(111) planes of diamond. 


Surface 

Number 

TersofT 

Brenner 

Plane 

of 

E b 

E b 


Adatoms 

(«V) 

(eV) 


i 

-7.4715 

-7.8397 

(100) 

2 

-7.5162 

-7.8413 


3 

-7.6803 

-7.9653 


1 

-3.5956 

-3.5210 

(111) 

2 

-4.8372 

-5.4031 

3 

-5.7485 

-5.9351 
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Figure 1. A top view for the (100) surface plane of diamond. Atomic sites in the 
top four layers (parallel to the exposed surface) are shown. Open large and open 
small circles represent atomic sites in the first and second layers, respectively. Solid 
large and solid small circles are atomic sites in the third and fourth layers. 
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Figure 2. A top view for the (111) surface plane of diamond. Atomic sites in the 
top four layers (parallel to the exposed surface) are shown. Open large and open 
small circles represent atomic sites in the first and second layers, respectively. Solid 
large and solid small circles are atomic sites in the third and fourth layers. In this 
case, atoms located in the second and third layers are superimposed. 
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Appendix 2. (2x1) Reconstructed Patterns of Diamond (100) Surface 


Simulation calculations were performed to investigate (2x1) reconstructed pat- 
terns for the (100) surface of diamond. Model potential functions developed recently 
for carbon species were used throughout this study. For a better analysis, two dif- 
ferent model functions, based on two and three-body interactions, were considered 
here. The first potential is the Tersoff function [1]. It has been shown that this 
function is able to reproduce correctly various bulk properties of diamond and of 
the graphitic plane. The second potential used in this study is the Brenner function 
[2] which is analytically similar to the Tersoff function, but parametrized differently. 
This function also, has been shown to produce acceptable results for properties of 
bulk diamond and the basal plane of graphite, as well as for some properties of small 
carbon clusters [2,7]. While the Tersoff function reproduces correctly the lattice con- 
stant of diamond, the Brenner potential at its minimum energy configuration gives 
a lattice constant value about 3% shorter. 

Relaxed and unrelaxed surface energies for low index planes of diamond for 
unreconstructed (lxl) patterns have recently been calculated using Tersoff and 
Brenner functions [8]. Unrelaxed surfaces created as abrupt terminations of the 
bulk have been found to undergo a multilayer relaxation process upon equilibra- 
tion. While a freshly created diamond (100) surface exhibits no LEED patterns, 
after heating to about 500 K in high vacuum it displays (lxl) patterns [3]. This 
surface, however, exhibits (2x1) patterns upon further heating to above 1300 K. 
Experimental evidence indicates that in this restructuring process the role of hy- 
drogen is significant. At the present time, it is generally believed that at high T, 
most of H is desorbed from the surface region, and the observed (2x1) LEED pat- 
terns are due to the formation of dimer pairs by analogy with the Si and Ge (100) 
surfaces. There are several theoretical calculations supporting the dimer formation 
model for the clean (100) surface of diamond. Recently, Mehandru and Anderson 
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[4] calculated the energy change associated with the dimerization process based 
on a six-layer-thick slab model. In their calculations they employed a tight bind- 
ing method (ASED-band) and demostrated that dimer formation is energetically 
favorable. Similarly, Bechstedt and Reichardt [5] using an energy minimization 
procedure based on a self-consistent tight-binding method, have concluded that a 
( 2 x 1 ) reconstructed ( 100 ) surface of diamond is energetically more stable than the 
unreconstructed (lxl) structure. Furthermore, considering a C 9 H 12 cluster (as a 
representative of the ( 100 ) surface) Verwoerd [ 6 ] has employed a MNDO procedure 
to calculate the minimum energy configuration of the top layer. His results also 
indicate that the dimer formation process is energetically favorable. 


In this investigation to determine relaxation energies as well as other struc- 
tural properties for differently restructured ( 2 x 1 ) patterns of the ( 100 ) surface of 
diamond, we carried out energy minimization calculations using both Tersoff and 
Brenner functions. A typical computational cell consisted of a 16-layer slab contain- 
ing about 256 atoms. Periodic boundary conditions were imposed on the system 
in two directions (i.e., x and y ) and the third direction, z, was left free to pro- 
vide exposed surfaces. Computational cell sizes in these calculations were chosen 
to satisfy minimum energy conditions for corresponding potential functions [ 1 , 2 ]. 
The relaxation energy, AE, was calculated as the difference in energies between a 
completly equilibrated system with a ( 2 x 1 ) structured surface and a system with 
an unrelaxed (lxl) face. Accordingly, the value of AE reflects the gain in stability 
as a result of restructuring and relaxation. 

The (100) surface of diamond with different (2x1) reconstruction patterns was 
taken into consideration. In addition to the (2x1) dimerized pattern, two different 
types of missing row models were also included in the calculations. In the first 
missing row (MR) structure, the top layer atoms located in every other row were 
removed from the system. A top view of the (MR) structure is shown in Figure 2 , 
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while Figure 1 depicts the unreconstructed (lxl) surface. The (MR) structure was 
obtained by removing alternating row atoms from an ideal (100) surface. (This sur- 
face structure is not identical with the structure obtained by removing alternating 
column atoms.) In this study, only the missing row model as shown in Figure 2 was 
investigated. In the second missing row model, those atoms removed to generate the 
(MR) structure, were placed on bridge sites along the rows formed by the remaining 
top layer atoms. Those displaced atoms now form a new layer which will be called 
the O’th layer. This model will be referred to as (MR+1) and its top view is depicted 
in Figure 3. A schematic top view of the (2x1) restructured dimer model is shown 
in Figure 4. This picture depicts dimers formed between the atoms within the same 
row. The other (2x1) dimer model which may be formed by pairing the column 
atoms, was found to be energetically less favorable. Therefore, in this investigation 
only the row-dimer model, as shown in Figure 4, was taken into, consideration. 


For the (2x1) dimer model, calculated relaxation energies along with the struc- 
tural parameters are compared in Table 1 with other available data from the lit- 
erature. While the gain in energy during the relaxation process calculated by the 
Tersoff function is smaller than the values of Mehandru and Anderson [4] and of 
Verwoerd [6], it is consistent with the value reported by Bechstedt and Reichardt [5]. 
On the other hand, the Brenner function produced a more negative AE value which 
shows, in a way, a better general agreement with reported values. The distance 
between two carbon atoms in a dimer, r dimer} calculated using the Tersoff function 
is in very good agreement with the value reported by Bechstedt and Reichardt [5]. 
On the other hand, the value of Tdimer obtained by the Brenner function is consis- 
tent with the reported value of Verwoerd [6]. Relatively speaking, both functions 
produced dimer distances in fair agreement with the literature values listed in Table 
1. In addition, both the Tersoff and Brenner functions predict symmetric dimers 
that agree with findings of Verwoerd [6] as well as with results of Mehandru and 
Anderson [4]. 
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Multilayer relaxation features for the ( 2 x 1 ) dimerized ( 100 ) surface calculated 
by the Tersoff and Brenner functions are qualitatively consistent with each other. 
SimilarOB to the (lxl) unstructured case reported earlier [8] for the dimerized 
(2x1) case also, interlayer spacings in the surface region exhibit alternating con- 
tractions and expansions with decreasing magnitudes as departing from the exposed 
surface. In this respect, results obtained here using the Tersoff as well as the Bren- 
ner functions are in general agreement with values of Mehandru and Anderson [ 4 ] . 
See Table 1. The change due to the relaxation in the interlayer spacing between 
the i and j planes is denoted by Ahjj and the corresponding percentage change 
with respect to the unrelaxed system is f>ij. The largest relaxation takes place in 
the first interlayer spacing which is the separation between the first and second 
surface planes parallel to the exposed surface. The value of Ahi2 obtained using 
the Tersoff function lies between the values reported by Mehandru and Anderson 
[ 4 ] and Verwoerd [6]. The Brenner function, on the other hand, produces a slightly 
larger value for Ahj2, but it is still in good agreement with the value of Verwoerd 
[6]. (See Table 1). Calculated values for AI123 and AI134 are very small. The Tersoff 
function predicts almost no relaxation for interlayer spacings between the third and 
fourth layers. When compared with reported values of Mehandru and Anderson [ 4 ], 
the Tersoff function, in general, predicts sligthly smaller values, while the Brenner 
function produces somewhat larger relaxation for the interlayer spacings. 

Calculated results obtained by the Tersoff and Brenner functions for the missing 
row models, (MR) and (MR+1), are presented in Table 2 . Tabulated values for AE 
indicate that for these missing row models the Brenner function provides more 
relaxations than the Tersoff potential. The same trend is also true for the dimer 
model. When those values are compared with relaxation energies given in Table 
1, the energetically least stable (2x1) restructured surface (for both functions) was 
found to be the (MR) model. The Brenner function predicts the dimer model 
to be the most favorable structure (among those three models considered in this 
investigation). It is about 0.18 eV per surface atom more stable than the (MR+ 1 ) 
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model. The Tersoff potential, on the other hand, predicts the (MR+1) structure to 
be the energetically most stable one. It was found to be 0.325 eV per surface atom 
more stable than the dimer model. Upon equilibration these missing row models 
also exhibit strong multilayer relaxation features. (See Table 2). For the (MR+1) 
model the relaxation patterns are qualitatively similar to the dimer model case 
with alternating contractions and expansions. For the (MR) model, however, the 
multilayer relaxation features are somewhat different. In this case, both first and 
second interlayer spacings shrink during the relaxation while the third interlayer 
spacing displays a relatively small amount of expansion. In addition to changes in 
interlayer spacings, relaxation in the case of these missing row models, also produces 
small lateral motions for some atoms. Those atoms located in the second layer 
were found to be displaced (with respect to their lattice positions) in a direction 
perpedicularly away from the missing row location. See Figures 2 and 3. For the 
(MR) case displacements were calculated as 0.07 and 0.11 A corresponding to the 
Tersoff and the Brenner functions; and for the (MR+1) case displacements were 
0.16 and 0.237 A, respectively. 

For the (100) surface of diamond both potential functions employed in this 
study predict that a (2x1) restructured dimer model is energetically more sta- 
ble than the (1x1) unrestructured surface. This result is consistent with earlier 
calculations found in the literature [4-6]. Among three different (2x1) patterns 
considered here, the Brenner function gives the dimer model as the energetically 
most favorable structure. The Tersoff function, on the other hand, predicts the 
(MS + 1) model to be the most stable configuration. So far, the dimer model has 
been proposed for the (2x1) reconstructed (100) surface of diamond by analogy 
with the Si(100) and Ge(100) surfaces. This model is very convincing and at the 
same time, it is well supported by theoretical calculations showing dimerization 
as an energetically favorable process. At the present time, however, there is not 
enough evidence indicating that the dimer model is the most favorable surface con- 
figuration leading to a (2x1) pattern. For that reason, perhaps, the choice of a 
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dimer model to explain the (2x1) LEED patterns may still need further substanti- 
ation. Despite many similarities, carbon in terms of structure, behaves somewhat 
differently than Si and Ge. For the bulk structure, for instance, the diamond phase 
is not the energetically most favorable one. More importantly, the (2x1) restruc- 
tured (100) surface of diamond does not exhibit the higher order reconstruction 
patterns which have been observered for Si and Ge cases [3,9]. In order to fully 
resolve this issue, in addition to new experiments, more accurate ab initio calcula- 
tions, perhaps for a number of different (2x1) patterns, are needed. Presently, we 
believe that the question “are there other (2x1) surface configurations with energies 
lower than the dimer model” remains unanswered for the (100) surface of diamond. 


Results obtained in this investigation support the stability of the dimer model 
over the (lxl) unreconstructed (100) surface of diamond. In fact, calculations 
indicate that all three (2x1) patterns which were considered in this work are ener- 
getically more stable than the unreconstructed (lxl) face. In general, calculated 
energies and structural properties for the dimer model were found to be in agreement 
with values reported in the literature. While the Brenner function, in accordance 
with general belief, predicts the dimer model to be the energetically most stable 
structure, the Tersoff function produced one of the missing row models, (MR+1), as 
the energetically most favorable form. This last result, however, should be weighed 
with caution. It may be an indication of a flaw in the potential function or we need 
to reconsider the assumption that the dimer model is the lowest energy configuration 
for the (2x1) restructured surface. 
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Table 1 . Structure and energy values for the (2x1) dimerized (100) surface of 
diamond (calculated by Tersoff and Brenner functions) compared with values from 
the literature. Relaxation energies, AE, are in eV per surface atom. The distance 
between two carbon atoms in a dimer is denoted by r<£j mer . The change in the inter- 
layer spacing between i and j layers is denoted by Ahij and 8{j is the correponding 
per cent change with respect to the unrelaxed structure. All distances are in A. 



Tersoff 

Brenner 

[4] 

[5] 

[6] 

AE 

-0.495 

-1.134 

-1.84 

(-1.26)* 

-0.41 

-1.55 

I* dimer 

1.542 

1.453 

1.58 

(1.40)* 

1.54 

1.434 

Ahj2 

-0.196 

-0.266 

-0.12 

(-0.14)* 


-0.24 

812 (%) 

-22.0 

-30.6 

-16.9 

(-21.4)* 



All23 

0.012 

0.124 

0.03 

(0.05)* 



*23 (%) 

1.34 

14.2 

6.7 

(9.0)* 



Ah 34 

0. 

-0.01 

-0.03 

(-0.03)* 



*34 (%) 

0. 

-1.0 

-3.4 

(-3.4)* 




* Spin paired results. See reference 4. 
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Table 2. Structure and energy values for (2x1) ’’Missing Row” models for the 
(100) surface of diamond calculated by Tersoff and Brenner functions. Descriptions 
of (MR) and (MR+1) models are given in the text. Relaxation energies, AE, are in 
eV per surface atom. The change in the interlayer spacing between i and j layers 
is denoted by Ahjj; and 6{j is the correponding per cent change with respect to the 
unrelaxed structure. The O’th layer is for the (MR+1) case, it corresponds to the 
new layer formed by the displaced first layer atoms (see the text). All distances are 

in A. 



Tersoff 

Brenner 

(MR) 

(MR+1) 

(MR) 

(MR+1) 

AE 

-0.391 

-0.820 

-0.704 

-0.955 

> 

o 


-0.157 


-0.252 

*01 (%) 


-17.63 


-29.05 

Ahi2 

-0.058 

0.183 

-0.098 

0.234 

*12 (%) 

-6.51 

20.52 

-11.29 

26.91 

All23 

-0.067 

-0.092 

-0.145 

-0.156 

*23 (%) 

-7.50 

-10.31 

-16.71 

-17.90 

Ah 34 

0.009 

0.010 

0.024 

0.019 

*34 (%) 

1.04 

1.22 

2.75 

2.23 
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Figure 1. A schematic top view for the (lxl) unreconstructed (100) surface of diamond. 

Positions of atoms located in the top three layers are shown. Large open 
circles represent top layer atoms. Small open circles indicate atoms in the 
second layer. Small solid circles are the third layer atoms. 
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Figure 2. Top view of a (2x1) reconstructed missing row model (MR) for the (100) 
surface of diamond. This structure was generated by removing carbon atoms 
located in alternating rows. Positions of atoms located in the top three layers 
are shown. Large open circles represent top layer atoms. Small open circles 
indicate atomic positions in the second layer and small solid circles are the 
third layer atoms. 
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Figure 3. Top view of a (2x1) reconstructed missing row model (MR+1). In this struc- 
ture, those atoms which were removed to generate the (MR) model, are in- 
troduced back into the system. They are placed on bridge sites of top layer 
atoms and shown by large open circles with a cross. These atoms constitute 
the O’th layer as indicated in the text. Large open circles represent atoms in 
the first layer. Small open circles indicate atoms in the second layer. Small 
solid circles are the third layer atoms. 
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Figure 4. Top view of the (2x1) restructured dimer model for the (100) surface of dia- 
mond. Positions of atoms located in the top three layers are shown. Large 
open circles represent top layer atoms. Small open circles indicate atoms in 
the second layer. Small solid circles are the third layer atoms. 
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Appendix 3. Crack Propagation and Tensile Behavior for 
the Basal Plane of Graphite 


The propagation of a surface crack was investigated at an atomistic level for 
a basal plane of graphite employing a computer simulation technique based on a 
molecular dynamics procedure. In addition to investigations of structural changes 
taking place during the crack propagation, in this study, the distribution of atomic 
stresses at the crack tip was also analyzed. Simulation calculations were carried out 
considering a potential function which has been derived recently by Tersoff [1] to 
represent energy- and structure-related properties of carbon systems. The Tersoff 
function, to date, has been used quite successfully in many studies. In addition to 
its ability to reproduce properties of a graphitic plane, this function can also predict 
many surface and bulk properties of the diamond phase. The total potential energy 
of a system of N particles in the Tersoff potential is expressed as a sum over atomic 
sites of the form: 

at 

t 

with 

1 N 

E < = 5 £ /.(rtfJIVMry) - lyVxfo)] (2) 

i (#0 

where, ry denotes the internuclear distance between particles i and j , and f c { r ij) 
represents the cut-off function which is given by: 

{ 1, if rij < R — D\ 

\ ~ ( T ij ~ D )l D \ » if R - D < ry < R + D] (3) 
0, if > R + D. 


The three-body part of the interactions is introduced via the function b{j , while Vr 
and Vr represent repulsive and attractive parts, respectively. These functions were 
defined by Tersoff as: 

Vn(rij) = A-exp(-A 1 r ij ) 
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VA{rij) = B-exp(-\ 2 rij) 

iy = (l +/3 n {,“)' 1/> “ 

N 

= ^ ] /c( r tfc)'5'(^ijfc)' ea: p[^3( 7 'ii ~ ^ifc) ] 

5(0iifc) = 1 + c 2 /d 2 — c 2 /[t£ 2 + (fe — cos^ijfc) 2 ] 

For carbon the Tersoff parameters are given as: >1 = 1393. 6(ey), B = 

346.74(eV), Aj = 3.4879(A), A 2 = 2.2119(A), /5 = 1.5724 x 10~ 7 , n = 0.72751, 
c = 38049., d = 4.3484, h = -0.57058, R = 1.95(A) and D = 0.15(A). 

Calculations were carried out considering a constant temperature molecular dynam- 
ics technique based on the Verlet algorithm [2]. Throughout this study a time step 
of 5. x 10 16 sec. and a reduced temperature of 0.25 were employed. Stresses for 
each atom in the system were evaluated as: 

_ 1 t dE < 

<7 “ /3 v dr} af3 

where, a a p denotes the atomic stress for the Cartesian components a and /3, and 
r) a fs represents the corresponding Lagrangian strain parameter [3]. The volume and 
the energy for the z’th atom are denoted by v and 25* , respectively. The model used 
in this study representing the basal plane of graphite contained 840 (30x28) atoms 
arranged in a rectangular two-dimensional array. First, a perfect lattice in a 2D 
honeycomb structure was brought to a static equilbrium. This system displays an 
almost elastic behavior and it resists any plastic deformation up to relatively high 
strain values. A crack was generated in this 2D system by removing 9 atoms from 
the surface region. A portion of this system with the crack is shown schematically 
in Figure 1. Then, a uniaxial load along the x direction was imposed on the system 
(which now bears a surface crack) in a stepwise fashion by progressively increasing 
the total length of the system in small increments. At any given time the actual 
elongation, e, is defined as: 

l-lo 
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where, l D and l denote the lengths of the system in x direction at the initial stage 
(i.e.j at time t = 0) and at a time t during the loading process. A periodic boundary 
condition was applied also in the x direction to provide continuity. In the perpen- 
dicular y direction, on the other hand, two exposed surfaces, one bearing the initial 
crack, were left intact. 

After the initial 1000 equilibration time steps, the system with a surface crack 
was elongated with small increments ( Ae = 0.005). Following each incremental elon- 
gation, the system was permitted to equilibrate for 50 time steps under the molec- 
ular dynamics code. In this second stage, the incremental loading and equilibration 
process was repeated 30 times corresponding to a total elongation of e = 0.161. 
During the final stage, no elongation was performed, but the system was simply 
permitted to relax for an extented period of time to observe the approach of the 
system to equilibrium. In Figure 2 a plot is presented displaying the elongation 
and the variation of the total energy of the system as a function of time steps. 
Except the very first few steps, the energy remains virtually constant during the 
first equilibration period. In the second stage, as expected, the energy increases 
progressively with the loading. Perhaps one may assume that at least in the begin- 
ning part of the second stage the system behaves elastically. As the existing crack 
propagates the stress is released and, as expected, the average energy of the system 
decreases. Figure 3 displays the variation of the average stress during the loading 
and equilibration stages. The phonon field generated in the lattice and its approach 
to equilibrium is visible in the final equilibration stage. These fluctuations in the 
stress value represent thermal effects [4]. Figures 4a through 4k display configura- 
tions of the system at various stages during the final relaxation period. As shown 
in these schematic drawings, the surface crack starts propagating after the 2500 th 
step. The solid circles in these figures represent 15 atoms with the largest stress 
values. In all cases, during the propagation period the tip of the crack exhibited 
itself as the high stress zone. This outcome is in general argeement with continuum 
studies. The crack propagates virtually perpedicular to the load direction. In this 
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study, no dislocation formation was observed in front of the crack tip. In order to 
investigate the size dependence of these results, calculations were repeated using 
a larger model containing 1680 (60x28) atoms. This larger system was found to 
produce results basically identical to those obtained by the smaller model. The 
time dependence of the propagation and the distribution of the high stress atoms 
as well as the crack tip geometry were all found to be very similar. 

In this case, only a very highly strained (e=0.161) system was taken into con- 
sideration. Despite the fact that present calculations correctly predict some of the 
contiuum results (such as the crack tip being the high stress zone), further calcula- 
tions with less strained systems are needed for a more complete analysis. Present 
calculations constitute the very first part of a general study on the simulation of 
crack propagation processes for covalently bonded systems. A proper represen- 
tation of atomic interactions in covalently bonded systems requires consideration 
of potential functions based on many-body interactions (like the Tersoff function 
employed here). One of the important features of the many-body interactions is 
to provide energetically stable configurations for open systems like the graphitic 
basal plane. In a crack propagation process, on the other hand, the role played by 
many-body forces constitutes a somewhat more complex academic problem; there- 
fore, it will be addressed separately. When these present results are compared with 
earlier calculations, the very short nature of the Tersoff function must be kept in 
mind. Furthermore, the ability of this function to simulate materials properties (in 
particular, those based on long range interactions) must be thoroughly analyzed. 
Another important point, in this respect, is related to dislocations. The absence of 
dislocations in this study may come from the brittle character of the model system 
considered here, or, more convincingly perhaps, it is due to the short range nature 
of the potential function employed. 
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Figure 1. A schematic representation for the initial configuration of the system with a 
surface crack. 
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